/** \file 
 \brief Accelerations by Verlet Cells.
 */

#ifndef YOCTO_CELL_DYN_INCLUDED
#define YOCTO_CELL_DYN_INCLUDED 1

#include "core/dynamics.hpp"

namespace FastMD  {
	
	//! Keep track of a particle
	struct CellNode
	{
		CellNode *next;
		size_t    index;
	};
	
	//! a Verlet Cell
	struct Cell
	{
		static const size_t MAX_LINKS = 3*3*3-1; //!< max #neighbors in 3D
		Real                V;                   //!< partial potential
		CellNode           *inside;              //!< first node in cell
		const size_t        links;               //!< #neighbors cells
		Cell               *link[MAX_LINKS];     //!< neighbors [0..links-1]
		const size_t        index;               //!< for info
		const U3D           id;                  //!< for info
		
		
		//! register a cell as a link, no multiple allowed
		/**
		 the cells are inserted by increasing address,
		 to avoid problems in the degenerate case of 2 cells
		 on one dimension.
		 */
		void subscribe( Cell *l ) throw();		
	private:
		Cell();~Cell();Cell( const Cell &); Cell & operator=( const Cell & );
	};
	
	//! Cells for one box
	class Cells 
	{
	public:
		
		//! a group of Cells
		struct Slice
		{
			size_t ini;
			size_t len;
			size_t end;
		};
		static const size_t MAX_THREADS = 128;
		
		static U3D CountFor( const V3D &box, const Real cutout) throw();
		
		Cells( const U3D &dim, const Frame &frame, const size_t num_threads );
		~Cells() throw();
		
		const U3D    size;  //!< #cells per dim
		const size_t count; //!< size.items()
		const U3D    side;  //!< size-1, for acceptable indices
		const size_t ystp;  //!< size.x
		const size_t zstp;  //!< size.x*size.y
		const V3D    ibox;  //!< frame box reciprocal
		
		inline Cell       & operator[]( const size_t index ) throw()       { assert(index<count); return cell[index]; }
		inline const Cell & operator[]( const size_t index ) const throw() { assert(index<count); return cell[index]; }
		
		void   dispatch( CellNode *node, const size_t nodes, const Vector &positions ) throw();
		
		inline U3D GetId( size_t index ) const throw()
		{
			U3D u;
			u.z = index / zstp; index -= u.z * zstp;
			u.y = index / ystp; index -= u.y * ystp;
			u.x = index;
			return u;
		}
		
	private:
		size_t       cmem;   //!< cmem >= count
		Cell        *cell;   //!< memory
	public:
		const size_t slices; //!< num slices
		Slice        slice[ MAX_THREADS ];
	private:
		YOCTO_DISABLE_COPY_AND_ASSIGN(Cells);
	};
	
	
	//! Common CellDynamics.
	class CellDynamicsBase : public Dynamics
	{
	public:
		virtual ~CellDynamicsBase() throw();
		const Real    cutout;
		void          build_cells( const Frame &frame );
		void          dispatch_only( const Vector &positions ) throw();
		const Cells & GetCells() const throw() {assert(cells!=NULL); return *cells; }
		
	private:
		size_t       nmem;  //!< nodes in memory
	protected:
		CellNode    *node;  //!< the nodes
		const size_t nodes; //!< #nodes
		Cells       *cells; //!< if allocated
	public:
		const        size_t num_threads; //!< required #threads
	protected:
		explicit CellDynamicsBase( const Frame &frame, const Real CutOut, const size_t nt);
		
		//! dispatch particles and set forces to 0
		void prolog( Frame &frame ) throw();
		
		//! update one slice of cells
		void work( size_t rank, Frame &frame , Potential &potential) throw();
		
		//! collect information and forces --> accelerations
		void epilog( Frame &frame, bool isothermal ) throw();
		
		
	private:
		uint64_t     wksp[ YOCTO_U64_FOR_ITEM(Cells) ];
		void release_nodes() throw();
		void release_cells() throw();
		YOCTO_DISABLE_COPY_AND_ASSIGN(CellDynamicsBase);
	};
	
	
	//! Single Threaded Dynamics.
	class CellDynamics : public CellDynamicsBase
	{
	public:
		explicit CellDynamics( const Frame &frame, const Real CutOut);
		virtual ~CellDynamics() throw();
		virtual void eval( Frame &frame , Potential &potential, bool isothermal) throw();
		
	private:
		YOCTO_DISABLE_COPY_AND_ASSIGN(CellDynamics);
		
	};
	
	
	
	
	
}

#endif
